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Abstract 

Background: Influenza A virus (lAV) infection-induced inflannnnatory regulatory networks (IRNs) are extremely 
complex and dynamic. Specific biological experiments for investigating the interactions between individual 
inflammatory factors cannot provide a detailed and insightful multidimensional view of IRNs. Recently, data from 
high-throughput technologies have permitted system-level analyses. The construction of large and cell-specific IRNs 
from high-throughput data is essential to understanding the pathogenesis of lAV infection. 

Results: In this study, we proposed a computational method, which combines nonlinear ordinary differential 
equation (ODE)-based optimization with mutual information, to construct a cell-specific optimized IRN during lAV 
infection by integrating gene expression data with a prior knowledge of network topology. Moreover, we used the 
average relative error and sensitivity analysis to evaluate the effectiveness of our proposed approach. Furthermore, 
from the optimized IRN, we confirmed 45 interactions between proteins in biological experiments and identified 37 
new regulatory interactions and 8 false positive interactions, including the following interactions: 111 (3 regulates 
TLR3, TLR3 regulates IFN-p and TNF regulates IL6. Most of these regulatory interactions are statistically significant by 
Z-statistic. The functional annotations of the optimized IRN demonstrated clearly that the defense response, 
immune response, response to wounding and regulation of cytokine production are the pivotal processes of 
lAV-induced inflammatory response. The pathway analysis results from the Kyoto Encyclopaedia of Genes and 
Genomes (KEGG) showed that 8 pathways are enriched significantly. The 5 pathways were validated by 
experiments, and 3 other pathways, including the intestinal immune network for IgA production, the cytosolic 
DNA-sensing pathway and the allograft rejection pathway, are the predicted novel pathways involved in the 
inflammatory response. 

Conclusions: Integration of knowledge-driven and data-driven methods allows us to construct an effective IRN 
during lAV infection. Based on the constructed network, we have identified new interactions among inflammatory 
factors and biological pathways. These findings provide new insight into our understanding of the molecular 
mechanisms in the inflammatory network in response to lAV infection. Further characterization and experimental 
validation of the interaction mechanisms identified from this study may lead to a novel therapeutic strategy for the 
control of infections and inflammatory responses. 
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Background 

Influenza A virus (lAV) infection is a worldwide public 
health threat [1,2]. lAV causes respiratory tract infec- 
tions and leads to inflammatory responses. Controlling 
the inflammatory response resulting from an lAV infec- 
tion is of great significance in reducing associated tissue 
damage. However, many biological experiments have 
demonstrated that lAV infection-induced inflammatory 
responses are extremely complicated and regulated by 
dynamic networks [3-5]. Specific biological experiments 
investigating the mechanisms of interactions among in- 
dividual inflammatory factors have not provided a suffi- 
ciently detailed and insightful multidimensional view of 
inflammatory regulatory networks (IRNs). We need to 
investigate the mechanisms at a system-level and from 
the network dynamics. Therefore, the construction of 
large and cell-specific inflammatory regulatory networks 
(IRNs) based on high-throughput data is essential for 
investigating the molecular mechanisms of inflamma- 
tory responses during lAV infection. 

Biological experiments have found that lAVs induce 
the expression of a number of inflammatory molecules 
and inflammatory cytokines and chemokines, such as 
IL27, IL32, IL6, TNF, IFNG, CXCLIO, CCL3, NOS2 and 
IL8 [6-9]. Furthermore, a number of studies have shown 
that the H5N1 viruses can induce increased gene tran- 
scription of pro-inflammatory cytokines, including 
CXCLIO, IFN-p, IL6, COX-2 (Cyclooxygenase-2) and 
CCL5 [9-12]. In particular, COX-2 is the primary medi- 
ator in protection against lAV infection [4] and has been 
shown to play a regulatory role in the induction of the 
H5Nl-mediated pro-inflammatory cascade [10,11]. It is 
important to further investigate the mechanisms of the 
inflammatory cascade downstream of COX-2 regulation 
that may be involved in H5N1 infection [13]. To our 
best knowledge, the studies on constructing a cell-specific 
IRN after lAV infection are limited, and an integrated and 
systematic analysis of the inflammatory cascade mediated 
by COX-2 that incorporates microarray data has not 
yet been reported. 

A number of different methods for inferring gene 
regulatory networks (GRNs) from high-throughput data 
have been proposed [14-20]. However, there are a few 
studies on the construction of dynamic signaling net- 
works based on stoichiometric approaches, discrete 
Boolean models, the fuzzy logic models, the integer pro- 
gramming method and the ordinary differential equa- 
tion (ODE)-based method [15,21-27]. No study has 
reported combining a prior knowledge of network topology 
with nonlinear optimization algorithms to identify the dy- 
namic regulatory network. In the process of reconstructing 
networks from expression data based on a priori knowledge 
of network topology, the most important steps are converting 
familiar network maps into mathematical models and fitting 



the available data into the networks structural para- 
meters. Recently, the rough topological structure of 
inflammatory networks with 2361 nodes and 63276 
edges in humans have been obtained, which provides 
a prelude to more detailed network analysis and 
mathematical modeling for an inflammatory network 
[28]. By combining information theory-based MI and 
nonlinear ODE-based optimization, in this study, we 
proposed a computational method to construct a cell- 
specific IRN mediated by COX-2 during lAV infec- 
tion. A differential evolution (DE) algorithm was used 
to optimize the network so that it best fits the experi- 
mental data. Furthermore, we performed a Kyoto En- 
cyclopaedia of Genes and Genomes (KEGG) pathway 
and gene ontology (GO) terms enrichment analysis on 
the optimized IRN to identify the underlying mecha- 
nisms during lAV infection. 

Methods 

The flowchart of our work is presented in Figure 1 
and mainly consists of six steps: constructing an initial 
IRN, simplifying the initial network, building a math- 
ematical model, estimating parameters in the model 
with the DE algorithm, significance test and sensitivity 
analysis for the regulations, and performing an enrich- 
ment analysis. 

Data collection and construction of the initial 
inflammatory regulatory network 

To construct a cell-specific IRN and investigate the 
mechanisms of the inflammatory cascade mediated by 
COX-2 in lAV infection, we selected 59 proteins, which 
are listed in Table 1, that are associated with the inflam- 
matory responses regulated by COX-2 based on the pub- 
lished literature [3,9,13,29,30]. The microarray data were 
retrieved from the Gene Expression Omnibus (GEO) 
database under the GEO accession number GSE28166 
[31,32]. This dataset contains 36 samples in total, with 3 
mock and infected replicates for each time point. In this 
study, the expression levels of complexes were the ave- 
rage of the gene expression levels of the members of the 
complex from the dataset. The expression level of lAV 
was obtained from the literature [32]. 

Network construction based on these 59 proteins was 
performed using Ingenuity Pathway Analysis (IPA) soft- 
ware (Ingenuity Systems, www.ingenuity.com ). The offi- 
cial symbol of each protein was imported into the IPA 
software. Through IPA analysis, we identified a total of 7 
networks based on functional connectivity. Three of 
these networks shared common proteins, and it was pos- 
sible to generate a merged network (data not shown). 
The merged network is very complicated and includes a 
few proteins that are not on our protein list. We pruned 
the network by removing the proteins that were not on 
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Constructing a initial inflammatory 
regulatory network based on priori 
knowledge 



Simplying network: removing the 
edges with independent correlation 
based on mutual information 



Building non-linear mathematical 
models for the networks 



Estimating parameters with DE 
algorithm 



NO ^^.^^ Stop ? 

Model fits experimental data 
well 

rYES 



Optimized network 



Significance test and sensitivity 
analysis for regulations 



Predicting regulatory interactions and 
preforming enrichment analysis 

Figure 1 The flowchart of our work. 



Simplifying the initial IRN based on conditional mutual 
information 

PCA-CMI, which was originally proposed by Zhang 
et al for inferring GRNs from gene expression data, 
considers the nonlinear dependence and topological 
structure of GRNs by employing a path consistency 
algorithm (PCA) based on conditional mutual infor- 
mation (CMI) [20]. In this study, we used the PCA- 
CMI method to distinguish direct (or causal) interac- 
tions from indirect associations. 

For a discrete variable X, the entropy H{X) is the 
measure of average uncertainty of variable X and can be 
defined by: 



(1) 



where p{x) is the probability of each discrete value x in X 
Mutual information (MI) measures the dependency 
between two variables (genes or proteins). For discrete 
variables X and Y, Ml is defined by the following equa- 
tion: 



(2) 



CMI measures conditional dependency between two 
variables given other variable(s). The CMI of variables X 
and y given Z is defined as: 

1{X,Y\Z) = - ^ ^(^,;,,,)log^^|MM_. (3) 



xeX,yGY,zeZ 



p{x\z)p{y\z) 



With the widely adopted Gaussian kernel probability 
density estimator, the equations (1), (2) and (3) can be 
easily calculated using the following equivalent equations 
[15,20]. 



H{X)=l\og{2ner\Cl 



/(x,r)=liogl^WNc(ni 

^ ' ' 2 ^ \C{X,Y)\ 



(4) 



(5) 



our list except for some common and important mole- 
cules, such as NFkB, IL12 (complex), p38 MAPK, JAK, 
STAT, IFN-a and IFN-|3. In addition, we integrated two 
molecules, lAV and COX-2, into the merged network. 
IL32, IL29, IL27, IL1|3 and IFN-a/|3/y have been 
reported to inhibit viral replication [3,33-37]. Therefore, 
we obtained our initial IRN comprising 51 proteins (or 
complexes) and 198 interactions. The network is 
depicted in Additional file 1. The full name of each pro- 
tein in the initial IRN is listed in Additional file 2. 



T(Y Y\7^-l ,^ JC{X,Z)\.\CiY,Z)\ 

^ ' '^"2 ^|c(z)|.|c(x, r,z)|' 



(6) 



where C is the covariance matrix of variable X, |C | is 
the determinant of matrix C, and n is the number of va- 
riables in C. 

A high MI value indicates that there is a close rela- 
tionship between the variables, while a low MI value 
implies variable independence. Similarly, a high CMI 
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Table 1 Total 59 proteins selected from the published 
literatures 



Table 1 Total 59 proteins selected from the published 
literatures (Continued) 



Gene Description 
symbol 



Type 



IL2RA Interleukin 2 receptor, alpha 



BPIFAl BPI fold containing family A, member 1 

CCL16 Chemokine ligand 16 

CCL2 Chemokine (C-C motif) ligand 2 

CCL22 Chemokine (C-C motif) ligand 22 

CCL5 Chemokine (C-C motif) ligand 5 

CCR9 Chemokine (C-C motif) receptor 9 

CD14 CD14 molecule 

CD27 CD27 molecule 

CMKLRl Chemokine-like receptor 1 

CSF2RA Colony stimulating factor 2 receptor, 
alpha 

CXCLIO Chemokine(C-X-C motif) ligand 10 

CXCL12 Chemokine(C-X-C motif) ligand 12 

CXCR3 Chemokine(C-X-C motif)receptor 3 

CXCR4 Chemokine(C-X-C motif)receptor 4 

CXCR5 Chemokine(C-X-C motif)receptor 5 

ERBB2 V-erb-b2 erythroblastic leukemia viral 
oncogene homolog 2 

FADD Fas (TNFRSF6)-associated via death 
domain 

FGF23 Fibroblast growth factor 23 

FGFRLl Fibroblast growth factor receptor-like 1 

HGF Hepatocyte growth factor (hepapoietin / 
scatter factor) 

HRG Histidine-rich glycoprotein 

IFNBl Interferon, beta 1, fibroblast 

IFNG Interferon, gamma 

IL12B Interleukin 12B 

IL15RA Interleukin 15 receptor, alpha 

111 6 Interleukin 16 

IL18BP Interleukin 18 binding protein 

ILip Interleukin 1, beta 

IL20RA Interleukin 20 receptor, alpha 

IL22RA2 Interleukin 22 receptor, alpha 2 

IL25 Interleukin 25 

IL27 Interleukin 27 

IL29 lnterleukin29(interferon, lambda 1) 



Other 

Cytokine 

Cytokine 

Cytokine 

Cytokine 

G-protein 
coupled receptor 

Transmembrane 
receptor 

Transmembrane 
receptor 

G-protein 
coupled receptor 

Transmembrane 
receptor 

Cytokine 

Cytokine 

G-protein 
coupled receptor 

G-protein 
coupled receptor 

G-protein 
coupled receptor 

Kinase 
Other 

Growth factor 

Transmembrane 
receptor 

Growth factor 

Other 

Cytokine 

Cytokine 

Cytokine 

Transmembrane 
receptor 

Cytokine 

Other 

Cytokine 

Transmembrane 
receptor 

Transmembrane 
receptor 

Cytokine 

Cytokine 

Other 



IL32 Interleukin 32 
IL6 Interleukin 6 (interferon, beta 2) 

IL6R Interleukin 6 receptor 

IL7 Interleukin 7 

KIT V-kit Hardy-Zuckerman 4 feline sarcoma 
viral oncogene homolog 

LTA Lymphotoxin alpha (TNF superfamily, 
member 1) 

LTBPl Latent transforming growth factor beta 
binding protein 1 

MET Met proto-oncogene 

MMP19 Matrix metallopeptidase 19 

NGF Nerve growth factor 

N0S2 Nitric oxide synthase 2, inducible 

NRG2 Neuregulin 2 

PECAMl Platelet/endothelial cell adhesion 
molecule 

SAAl Serum amyloid Al 

SELL Selectin L 

SMAD9 SMAD family member 9 

TGFBl Transforming growth factor, beta 1 
TIMP4 TIMP metallopeptidase inhibitor 4 
TLR3 Toll-like receptor 3 

TMEFF2 Transmembrane with EGF-like and two 
follistatin-like domains 2 

TNF Tumor necrosis factor 

TNFSFIO Tumor necrosis factor (ligand) superfamily, 
member 10 

TNFSF12 Tumor necrosis factor (ligand) superfamily, 
member 12 

TNFSF14 Tumor necrosis factor (ligand) superfamily, 
member 14 

The gene symbol, description and type corresponding to each protein 
are given. 



indicates that there is a close relationship between the 
variables X and Y given variable Z, while a low CMI 
value represents independence between genes. If the MI 
or CMI value is smaller than a given threshold 6, the 
edge between the two proteins is deleted for the inde- 
pendence (See the detailed procedure of PCA-CMI in 
[20]). 

We simplified the initial IRN based on PCA-CML We 
deleted the edges of the initial IRN with independent 



Transmembrane 
receptor 

Other 

Cytokine 

Transmembrane 
receptor 

Cytokine 

Kinase 

Cytokine 

Other 

Kinase 
Peptidase 
Growth factor 
Enzyme 
Growth factor 
Other 

Transporter 
Other 

Transcription 
regulator 

Growth factor 

Other 

Transmembrane 
receptor 

Other 



Cytokine 
Cytokine 

Cytokine 

Cytokine 
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correlations recursively (i.e., from low to high order of 
independent correlation until no edge can be deleted). 
The threshold value 6 of MI and CMI was set to 0.1, 
which was statistically tested by Z-statistic [15,20,38] 
(Figure 2). The simplified network, which contains 50 
nodes and 142 directed edges, is depicted in Figure 3. 

Mathematical model of the network 

To further obtain the more simplified IRN, we built 
nonlinear ordinary differential equations (ODEs) to 
model the reaction kinetics of the regulatory network. 
The ODEs describe the relationship between the reac- 
tion rate and the concentrations of the reactants. The 
change in concentration of a reactant is characterized 
by a function that takes the regulatory influence (acti- 
vation or inhibition) of other reactants into account. 
The general form of nonlinear ODEs is described as 
follows: 



dt 



where Xi is the concentration of species /, fi is a nonlinear 
function, m is the number of species in the system, kij is 
the kinetic parameter with ;e{l, 2, ... , m} and di is the 
degradation rate of species /. 

Based on the law of mass action and Hill functions, 
the nonlinear ODEs including 50 equations and 192 kin- 
etic parameters were built. All equations and their expla- 
nations and the initial concentrations of proteins are 
listed in Additional file 3. 



0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 

Z-statistic 

Figure 2 The significance test for the selected threshold value 
of the Ml and CMI. The x-axis is the Z-statistic values of the deleted 
edges when the threshold value of Ml and CMI is set to be 0.1 . The 
y-axis represents the number of edges whose Z-statistic fall into the 
corresponding bins. The blue dashed line is the inverse cumulative 
distribution function of A/(0,1) when using significant level a=0.1. 



Estimation of the kinetic parameters in the model with 
the DE algorithm 

The parameters in our ODEs can be classified into 
two categories of regulatory parameters: parameters 
representing activation or inhibition relations and deg- 
radation parameters representing the degradation of 
individual biomolecular species. 

The problem that identifies the kinetic parameters in 
the model can be converted into the following nonlinear 
optimization problem, which is the minimization of the 
error between the simulation values in our model and 
the experimental data. 

N M 



i=l j=l 

where xf^'^itj) and xf^{tp K) are the experiment data and 
simulation value of the species / at tj time point respec- 
tively. o)i = l/ max Xi^^'^iytj))^ , K is the parameter set 

consisting of all the parameters in the model, N is the 
number of species and M is the number of time points 
in the biological experiments. 

A wide variety of global optimization techniques have 
been developed to address nonlinear optimization prob- 
lems [39-42]. The DE algorithm, proposed by Storn and 
Price [42], is a very successful and powerful population- 
based stochastic search technique for solving global 
optimization problems and has been widely applied in 
many scientific and engineering fields [43,44]. The DE 
algorithm is described as follows: 

Step 1. Initialization: Generate random initial popula- 
tion Xg={ Xi^G > ^ZG Xm,g }> where X^g={ -^/i.G > -^/2. 
G XiD,Q }, N is the population size , G is the gene- 
ration, X/ G is a rate constant set and D is the amount of 
the kinetic parameters. 

Step 2. Genetic operation: 

1) Mutation: Vi^G+i=Xri,G-^F{XrZG- ^rs^c)^ ri, r2> 
r3e{l,...,7V}\{i},FG[0,2]. 

2) Crossover: Ui^G+i=(^uG+i^"MDi,G+i) 



U 



f Vji^G+i. if {rand{b(j))<CR) or ; = rnbr{i) 
I Xy/^G+i^if {rand{b(j)) > CR) ov j^rnbr{i) ' 

1,2, "'D, 



where rand{b{j)) is the ;th evaluation of a uniform 
random number in [0,1], CR is the crossover con- 
stant in [0,1] and rnbr{i) is a random indexes in 

Step 3. Selection: If /(/7,,g) <M,g), then X,,g-h i = ^/,g, 
else X^,G+ 1 =^/,G> where /is the objection function. 



Jin and Zou BMC Systems Biology 2013, 7:105 
httpy/www.biomedcentral.com/l 752-0509/7/1 05 



Page 6 of 14 




Jin and Zou BMC Systems Biology 2013, 7:105 
httpy/www.biomedcentral.com/l 752-0509/7/1 05 



Page 7 of 14 



Average relative error 

The average relative error (ARE) is defined as follows. 



ARE 



i^|r--(^.)-rexp(^.)| 

1=1 



where and Y^'^'^it^ are the simulation and experi- 

ment values of the protein at time point ti and n is the 
number of samples. In this study, n^6. 

Sensitivity analysis 

Sensitivity analysis is a useful way to investigate the ef- 
fects of parameters variations on changes in the model 
outputs. We formulate the sensitivity Si{t) of parameter 
P at time t as follows: 

'^^'> - 0,{t) ' P{t) 

«(|Oi(P + AP, t)-Oi{P, t)\/Ot{P, t))l{AP/P), 

where 0/(t) is the /-th model output at time t, P is the 
parameter, AP is a small perturbation of P, 

Then we define the sensitivity Si of the i-th model out- 
put with respect to parameter P blow. 



where n is the number of samples (time points). 
Enrichment analysis 

We conducted a functional enrichment analysis for the 
network based on GO Biological Processes (BP) terms 
and the KEGG pathway with the DAVID bioinformatics 
database [45]. The enrichment significance was deter- 
mined by the DAVID tool. The P-values were then 
corrected for the false discovery rate (FDR). In this 
study, all the proteins other than lAV in the network are 
mapped with the DAVID database. For the complex, one 
member of the complex was mapped. The criterion for 
statistically significant enrichment was an FDR adjusted 
p-value less than 0.002. 

Results 

The optimized IRN based on the experimental data 

The initial and simplified IRNs (Figure 3) were 
constructed using IPA software and the PCA-CMI algo- 
rithm (See the Methods), respectively. To further 
optimize the network according to the experimental data, 
we first estimated all parameters in our nonlinear ODEs 
by the DE algorithm (See the Methods). The DE algo- 
rithm was carried out ten times, and the best parameter 
set was obtained, which is listed at Additional file 4: 
Table S2. 



Second, we further deleted some nodes and edges to 
simplify the IRN according to the following rules. If the 
optimal value of the kinetic parameter kij was zero, we 
deleted the directed edge, which indicates that biomole- 
cular ; does not regulate biomolecular / in the network. 
Furthermore, if there was no edge to connect with biomo- 
lecular /, we deleted the node / in the network. Finally, if 
the node / has been deleted in the network, the degra- 
dation rate di was set to zero in the numerical simulation. 
The optimized IRN is shown in Figure 4. 

Based on the optimal parameters, we performed a nu- 
merical simulation for all nodes in the network for com- 
parison with the experimental data. The dynamical 
processes of 8 key proteins are plotted in Figure 5 and 
those of other proteins are displayed in Additional file 5. 
The average relative errors (AREs) of the 98% proteins are 
less than 0.3, and those of the 2% proteins are within the 
interval [0.3, 0.7] (Figure 6). These results indicated the fi- 
delity of the obtained IRN. In addition, from the dynam- 
ical viewpoint, sensitivity analysis of the ODE models is 
very important to quantify the reliability of the parameters 
(regulatory strength between two genes) in the model (See 
the Methods). The results of the sensitivity analysis 
showed that the concentrations of the proteins are not 
sensitive to the perturbation of parameters (Figure 7), 
which indicating the reliability of the obtained IRN. 

Prediction of regulatory interactions in IRN 

Among the regulatory interactions in the optimized net- 
work, 45 interactions have been reported in the literature 
and are represented by red lines in Figure 4. In addition, 37 
new regulatory interactions have been predicted from the 
network and are denoted by black lines in Figure 4. Fur- 
thermore, the statistical significance of these regulations 
between paired proteins was tested using the method 
presented in the literature [15,20]. The significant and 
non-significant regulations were denoted by thick and thin 
lines in Figure 4, respectively. The number of significant 
and non-significant regulations was summarized in Table 2. 
The results demonstrated that most of the predicted regu- 
latory interactions, which are the same as the validated 
experimental interactions, are statistically significant. 

The presence of false positive interactions is a common 
problem in inferring a network. One source of false posi- 
tive interactions is indirect effects (i.e., in a cascade A^ 
B^ C and A^ C, protein A regulates C and may be me- 
diated by B, so the direct regulatory interaction A^ C 
may be a false positive interaction). Comparing the opti- 
mized IRN with the initial IRN, we have also identified 8 
false positive interactions, which are shown by dashed 
lines in Figure 3. For example, the interactions involving 
the regulation of IFN-|3 by TLR3 and IL6 by TNF take ef- 
fects during lAV infection through other chemical mole- 
cules. In our work, we have found that TLR3 regulates 
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Figure 4 The optimized inflammatory regulatory networlc. The lines ended with arrows and bars denote positive and negative regulatory 
interactions, respectively.The lines without arrows or bars represent binding interactions. The dashed lines indicated the false positive interactions 
identified by the proposed method. The red and black lines stand for the regulatory interactions which have been validated by biological 
experiments and are predicted by the proposed method, respectively. The significant and non-significant regulations are denoted by thick and 
thin lines, respectively. 



IFN-p through NFkB signaling, which is consistent with 
previous findings. The TLR3-induced NFkB signaling 
pathway is triggered by the virus, and NFkB regulates 
expression of the proinflammatory molecule IFN-p in the 
immune responses [46]. We have also found that TNF reg- 
ulates IL6 mediated by the activation of CCL2 or CD 14. 
The interactions in the optimized network are further 
classified, and detailed descriptions are presented in 
Additional file 6: Table S3. 

Identification of the important biological processes and 
pathways 

To gain further insight into the biological interpretation 
of the optimized IRN during lAV infection, we have 



performed Biological Process (BP) terms and a KEGG 
pathway enrichment analysis of all the species other than 
lAV in the optimized IRN (See the Methods). Functional 
enrichment analysis of the species was conducted using 
DAVID [45]. The annotation analysis shows enrichment 
in BP involved in the defense response, inflammatory re- 
sponse, immune response and regulation of cytokine 
production. The top 10 enriched BP are listed in Table 3, 
and all detailed lists of the significantly enriched BP 
terms (FDR<0.002) are available in Additional file 7. In 
addition, network ontology analysis (NO A) which anno- 
tates biological networks [47], was used to analyze the 
enriched functions of the optimized IRN. The results of 
NO A are shown in Additional file 7. The functions of the 
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Figure 5 Comparisons between tiie numerical simulation results and experimental data of lAV, IL32, IFN-p, TLR3, CCL16, CD40, TNF, NFkB. 

The blue and red lines denote the experiment and simulation results, respectively. The stars represent experimental data at each time point. The 
experimental errors are also plotted as short bars at each time point. 
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ARE 

Figure 6 The distribution of the average relative errors (ARBs) 
for the numerical simulation of the proteins in the optimized 
network. The y-axis represents the number of proteins whose AREs 
fall into the corresponding bins. 
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Figure 7 The sensitivity analysis for the parameters in the ODE 
models. The x-axis is the outputs (proteins) in the model and the 
y-axis is the calculated sensitivity. 
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Table 2 The number of significant and non-significant 
regulations 





Validated 


Predicted 


Significant 


24 


21 


Non-significant 


21 


16 



The second and third columns are the number of regulations validated by 
biological experiments and predicted by our proposed method, respectively. 



optimized lAV-induced IRN are enriched in the response 
to stimulus, immune system process, inflammatory re- 
sponse, response to wounding and positive regulation of 
cellular processes, which are similar to the results by 
employing DAVID. The functional annotations of the op- 
timized IRN reveal that the reconstructed IRN functions 
reasonably well and they reflect the defense response, 
immune response, response to wounding and regulation 
of cytokine production are the important processes of 
lAV-induced inflammatory response. 

The KEGG pathway enrichment analysis have identified 
a total of 8 pathways significantly enriched in this network 
(FDR<0.002). These pathways are shown in Table 4. 
Among them, the Toll-like receptors, the RIG-I-like recep- 
tors and the NOD-like receptors have been shown to be 
the main pattern-recognition receptors (PRR) by which 
the innate immune system recognizes the influenza virus 
infection [48]. Moreover, the NOD-like receptors play a 
primary role in host defence against invading pathogens 
and regulating NFkB signalling, ILlp production, and cell 
death, indicating that they are crucial to the pathogenesis 
of a variety of inflammatory human diseases [49]. The 



Table 3 Top 10 significantly enriched GO terms 



GO term (Biological Processes) 


Count 


FDR 


GO:0006952~defense response 


24 


7.01 e- 
20 


GO:0006954~inflammatory response 


18 


9.38e- 
16 


GO:0009617~response to bacterium 


15 


2.53e- 
14 


GO:0006955~immune response 


21 


2.65e- 
14 


GO:000961 1 -response to wounding 


18 


3.56e- 
12 


GO:0006935~chemotaxis 


12 


2.35e- 
10 


GO:0042330~taxis 


12 


2.35e- 
10 


GO:0002684~positive regulation of immune system 
process 


13 


6.24e- 
10 


GO:0051240~positive regulation of multicellular 
organismal process 


13 


8.39e- 
10 


GO:0001817~regulation of cytokine production 


12 


9.18e- 
10 



The second column is the number of proteins from our protein list in the 
corresponding GO term. 



Table 4 Significantly enriched KEGG pathways 



KEGG pathway term 


Count 


FDR 


References 


H sa 04060:Cyto ki n e-cyto ki n e rece pto r 
interaction 


23 


z.4ye- 1 y 


[49,50] 


Hsa04620:Toll-like receptor signaling 
pathway 


14 


5.1 1 e-1 2 


f/i "7 cn col 

L4/,bU,bzJ 


Hsa04672:lntestinal immune network for 
IgA production 


8 


1 .07e-05 


Predicted 


Hsa04630:Jak-STAT signaling pathway 


11 


1 .30e-05 


[36,49,50] 


Hsa04622:RIG-l-like receptor signaling 
pathway 


8 


1.51e-04 


[47,53-55] 


Hsa04623:Cytosolic DNA-sensing 
pathway 


7 


7.10e-04 


Predicted 


Hsa04621:NOD-like receptor signaling 
pathway 


7 


0.001461 


[47,48,50] 


Hsa05330:Allograft rejection 


6 


0.001944 


Predicted 



The second column is the number of proteins from our protein list in the 
corresponding GO term. In the last column, the corresponding references of 
the pathways involved in influenza A virus infection are given, and the terms 
"predicted" indicate that the pathways have not been reported by 
previous literatures. 



cytokine-cytokine receptor interaction and Jak-STAT sig- 
nalling pathway are also well known antiviral response 
pathways [50,51]. 

Three additional identified pathways have not been 
demonstrated to be associated with lAV infection. The 
intestinal immune network for IgA production signifi- 
cantly enriched (FDR=1.07e-05). Some researchers have 
reported that serum IgA is an inflammatory antibody 
that interacts with FcaRI on effector immune cells and 
may function as a second line of defence by eliminating 
pathogens that have breached the mucosal surface 
[56,57]. The detection of cytosolic DNA is related to the 
induction of IFN-a/|3 and other pro-inflammatory cyto- 
kines [58-61]. Cytosolic DNA has also been shown to ac- 
tivate the TBKl, IRF3 and the caspase-l-dependent 
maturation of IL-1|3 and IL-18 [58,62]. Allograft rejec- 
tion is also enriched significantly. Some authors have 
reported that influenza infections are associated with 
allograft rejection, but there is no evidence that lAV trig- 
gers the acute allograft rejection episodes [63-65]. In our 
results, under the stimulation of lAV, the allograft rejec- 
tion pathway is significantly enriched. These three path- 
ways lack literature support but may be promising novel 
pathways and need the experimental validation. 

Discussion 

The induction of pro-inflammatory cytokines such as 
COX-2, TNF, IFNs, IL27 and CXCLIO is essential for 
the host immune response during virus infection, but 
inappropriately sustained induction causes cytokine- 
storms, which are associated with a wide variety of infec- 
tious diseases [66,67]. Because of the complexity of the 
inflammatory response, it is necessary to study the 
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underlying mechanisms of inflammatory response based 
on a network approach. In this study, we proposed a 
nonlinear ODE-model based computational method to 
construct a cell-speciflc IRN during lAV infection. The 
main contributions of this study include three aspects. 
First, we built the large-scaled nonlinear ODE model of 
the network including 50 equations and 192 kinetic pa- 
rameters. Most of model-based studies for inferring net- 
works are based on linear ODE models or discrete 
models [14-16], and these linear ODEs are approximated 
by difference equations or the steady-state assumption, 
which are easily solved by classical optimization algo- 
rithms or software. However, the regulatory interactions 
in real biological networks are often non-linear. There- 
fore, the non-linear ODE model can better describe the 
complicated regulatory networks. The comparison study 
for the advantage of involving nonlinear items in the 
model was also performed by using linear ODE model 
to describe the regulatory network. The AREs in the 
linear model exhibited significantly higher values than 
those in the nonlinear model (Additional file 8: Figure SI, 
P-value<0.001). These results indicated that the non-linear 
ODE model can better describe the complicated regulatory 
networks. Second, we combined the DE algorithm with 
a priori knowledge to refine the nonlinear ODEs and 
solve the nonlinear optimization problem derived from 
constructing the network. This nonlinear optimization 
problem is difficult to solve using classical optimization 
algorithms because of high nonlinearity and no explicit 
expression. Although DE algorithm is a published stochas- 
tic search technique, it is a repeated process from the 
model to optimization and then from improved model to 
optimization. If the model is not correct, the best 
optimization algorithm is also useless. Our nonlinear ODE 
model has been repeatedly adjusted. Finally, global errors 
that reflect the effectiveness of fitting the reconstructed 
network to experimental data are presented. In most stu- 
dies based on the linear model systems, they did not pro- 
vide the errors or only gave the residual errors (local 
errors) that cannot quantify the real error between the 
networks and the experimental data. 

Because our proposed method integrated gene expres- 
sion data with a priori knowledge of topological struc- 
ture from literature and IPA software, it cannot compare 
with the published purely data-driven methods to evalu- 
ate the predictive results. However, these published ex- 
cellent works may help us to find a more appropriate 
way to evaluate the approaches that combined the ex- 
perimental data and a priori knowledge in the future. 

An increasing number of researchers have focused on 
the gene expression profile of host cells infected by in- 
fluenza virus [68-70]. However, most reports involve a 
single gene or pathway [52,53,71]. Few studies have fo- 
cused on the systematic analysis of the regulation of the 



cell-signaling cascade by lAV. To understand the global 
regulatory mechanisms of the inflammatory response 
during lAV infection, we conducted a pathway en- 
richment analysis of the optimal IRN with the KEGG 
database. From our results, a few host cellular signaling 
pathways stimulated by lAV infection have been identified. 
Some of these signaling pathways are critical to the innate 
immune response of the host cell against influenza virus, 
such as the Toll-like receptor, the RIG-I-like receptor and 
the NOD-like receptor pathways [48,54]. The activation of 
the TLR signaling pathway results in the stimulation of 
both innate and adaptive immune responses, and TLR 
agonists may represent an effective and broad-spectrum 
antiviral strategy to combat influenza viruses [71]. Several 
virus-encoded components that antagonize RLR signalling 
interact with and inhibit the IFN-a/|3 activation pathway 
using both RNA-dependent and RNA-independent me- 
chanisms [55,72]. 

Among the three novel pathways identified in our study, 
the functions of IgA have been studied [56,57]. Secretory 
immunoglobulin A (SIgA) is the major antibody isotype 
present in mucosal secretions and has many fiinctional at- 
tributes, both direct and indirect, serving to prevent in- 
fective agents such as bacteria and viruses from breaching 
the mucosal barrier [42]. Many DNA-sensors have been 
reported, such as IFI16, RNA Polymerase III, DAI, AIM2, 
NLRP3, LRRFIPl and DDX9/36. They play an important 
role in IFN-a/p and cytokine production [54,58,73]. For 
example, IFI16 can induce the inflammasome in response 
to Kaposi's sarcoma-associated herpesvirus infection and act 
as a mediator of the anti-inflammatory actions of type I IFNs 
[73-76]. AIM2 triggers the assembly of the inflammasome, 
culminating in caspase-1 activation, IL-1|3 maturation and 
pyroptotic cell death [77]. LRRFIPl has been shown to 
contribute to the production of IFN-p induced by vesicu- 
lar stomatitis virus (VSV) and Listeria monocytogenes in 
macrophages [78]. This evidence indicates that DNA 
sensors play an important role in virus infection. 
However, these results need further biological experimental 
verification. 

It should be noted that the protein activity profiles are 
substituted with the corresponding gene expression levels 
in the computation in the study because the protein ac- 
tivity profile data have not been easily obtained thus far. 
Therefore, there may be some discrepancy when mo- 
delling the network. In addition, the network we 
constructed does not involve RNA components, such as 
target mRNAs, micro-RNAs (mi-RNAs) or other non- 
coding RNAs, which may also modulate signals at many 
steps. Recent studies have provided evidence of a poten- 
tial role for long non-coding RNAs (IncRNAs) in regula- 
ting inflammatory gene expression [79,80]. Emerging 
evidence shows that mi-RNAs have been clearly impli- 
cated in the regulation of inflammatory responses 
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[81-85]. To better understand the molecular mechanisms 
of the inflammatory response during lAV infection, it 
requires the challenging process of constructing inflam- 
matory regulatory networks by integrating different types 
of data, such as gene expression data, protein activity 
profiles, mi-RNAs expression profiles and Chip-seq data. 

Conclusions 

A cell-specific IRN in lAV infection has been cons- 
tructed based on the proposed method. Furthermore, 
37 new regulatory interactions were predicted and 8 
false positive interactions of IRN and 3 novel pathways 
were identified in the study. These new findings pro- 
vide insight into our understanding of the mechanism 
of inflammatory response in lAV infection. Under- 
standing the pivotal role of signaling pathways during 
lAV infection may provide new insight into therapeutic 
strategies for the control of virus infection and inflam- 
matory response. Our findings also have significant 
implications on the development of biomarkers for in- 
fectious disease. 
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